Impurity model for non-equilibrium steady states 
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We propose an out-of-equilibrium impurity model for the dynamical mean-field description of the 
Hubbard model driven by a finite electric field. The out-of-equilibrium impurity environment is 
decomposed into a collection of equilibrium reservoirs at different chemical potentials. We discuss 
the validity of the impurity model and propose a non perturbative method, based on a quantum 
Monte Carlo solver, which provides the steady-state solutions of the impurity and original lattice 
problems. We discuss the relevance of this approach for other non-equilibrium steady-state contexts. 



The understanding and techniques of condensed mat- 
ter theory are such that it is possible nowadays to make 
reliable quantitative predictions on equilibrium prop- 
erties of many materials, including strongly-correlated 
many-body systems. Despite those many progress, very 
little is know in comparison for correlated systems far 
from equilibrium. Even for simple measurements such 
as current-voltage characteristics, one can only hope 
for some qualitative theoretical understanding. In this 
perspective, the study of non-equilibrium steady states 
(NESS) is paramount for understanding which of the 
equilibrium concepts, techniques and results can be 
adapted to non-equilibrium situations. 

The adaptation of equilibrium numerical methods 
to time-dependent situations {e.g. exact diagonaliza- 
tion, density matrix renormalization group, diagram- 
matic Monte Carlo, dynamical mean-field theory, vari- 
ational methods) enable the study of complex transients 
such as the dynamics after a quench in temperature or in- 
teraction. However, the limitations in time of such time- 
dependent methods (e.g. size or Poincare recurrence ef- 
fects for close systems, growth of entanglement entropy, 
sign problem) usually make them unsuitable for the study 
of steady-state physics. 

Recent theoretical efforts have been done to address 
the non-equilibrium steady-state physics of quantum dots 
driven by a constant voltage by developing steady-state 
formalisms [l|, 0] that bypass the transient dynamics. 
Amongst several promising steady-state techniques Q, 
a remarkable leap forward has been recently achieved 
by Han and Heary who developed a non-equilibrium 
steady-state impurity solver based on a Matsubara-like 
formalism and a Hirsch-Fye algorithm (NESS-HF) [3,@- 
Progress has also been made for describing the non- 
equilibrium steady-state dynamics of correlated electrons 
on finite dimensional lattices [(J 0] ■ In the context of cor- 
related electrons driven by a constant electric field, we 
recently proposed a dynamical mean-field theory solu- 
tion of the non-equilibrium steady state (NESS-DMFT) 
based on a gauge- invariant Keldysh formalism [1, Q . 

In this work, we give substance to the NESS physics 
of such lattice problems by focusing on their local mean- 



field description (the so called impurity problem), estab- 
lishing the connection with the out-of-equilibrium physics 
of quantum dots. In the framework of the electric-field- 
driven Hubbard model, we model the environment of the 
related impurity problem by a collection of equilibrium 
leads at different chemical potentials. We show how non- 
equilibrium properties on the lattice side, such as the 
energy distribution function or the dissipation, translate 
on the impurity side. Simplifying the impurity model by 
truncating it to the first relevant leads, we solve it numer- 
ically via a generalized NESS-HF solver a la Han-Heary, 
and we obtain the corresponding NESS-DMFT solution 
of the electric-field-drivcn Hubbard model. We compare 
the results with the ones obtained by solving the impurity 
to the second order in the interaction with the so-called 
iterated perturbation theory (IPT). 

Lattice model. We consider the Hubbard model on a 
d = 2 square lattice. It is driven out of equilibrium by 
a static and uniform electric field set along the x-axis of 
the lattice: E = Eu x with E > 0. The corresponding 
Lagangian reads (we set H = 1 and use the conventions 

Cs = ^2 Cia [idt - <&(*)] C ia ~U^2 CifCifCilCil 
+ Cirhj^^Cjo + conj. 

where Ci a and ci a are the Grassmann fields representing 
an electron at site i with spin u G {til}- The Hubbard 
Z7 term accounts for the on-site Coulombic interaction 
and Uj = (a/27r) 2 Jdk e lk ' XiJ e(k) is the hopping ampli- 
tude between two nearest neighbors distant of a: e(k) = 
eo[cos(k x a) + cos(k y a)]. ctij(t) = qf *dx ■ A(t,x) are 
the Peierls phase factors, q is the charge of the elec- 
trons and A ((f>) is the vector (scalar) potential: E = 
-V^-ftA. 

To allow for non-trivial steady states, we couple 
the system to a thermostat composed of independent 
reservoirs in equilibrium at temperature T: C s b — 
7j2io-i e 10i( - t * > biaiCia+ conj. . The b^i's are non-interacting 
electrons, and I labels their energy level. The phases 
6i(t) = J dt' 4>i(t') ensure the U(l) gauge invariance and 
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we work with a particle-hole symmetric system by sim- 
ply considering half-filled reservoirs with a flat density 
of states of bandwidth W. The thermostat introduces a 
local (gauge-invariant) retarded hybridization — ir where 
T = 7 2 /W sets the dissipation rate. The corresponding 
Keldysh component of the hybridization is fixed by the 
fluctuation-dissipation theorem. We take the tempera- 
ture T of the thermostat to be the lowest energy scale 
(fee = 1). We restrict ourselves to the paramagnetic so- 
lution and drop the spin indices. In the numerics, we use 
q = a = 1 and measure energies in units of eo- 

Dynamical mean-field approach. In equilibrium, the 
self-energy of a <i-dimensional Hubbard model reduces 
to the one of a dimensionless Anderson Impurity Model 
in the limit of infinite connectivity. Dynamical mean- 
field theory (DMFT) uses this result to compute non- 
perturbative mean-field solutions of finite dimensional 
lattice problems by solving auxiliary impurity problems 
defined self-consistently |l0|. Here, we propose to gener- 
alize this mapping for non-equilibrium steady states. 

The Keldysh action of a single-band Hubbard- U im- 
purity in a generic steady state reads 



lattice 



ab 



s = EE Utdt'<z(t)g^ ab (t-t')cl(t') 



J> At Uc a t (t)c a t (t)cl(t)cl(t) , 



(2) 



where the indices a, b = ± refer to the forward and back- 
ward branch of the Keldysh contour. The Q§ b are the im- 
purity non-interacting Green's functions (often referred 
as the Weiss effective fields). The two quantities that 
fully describe the non-equilibrium steady-state physics 
of the impurity are its density of states p(e) and its en- 
ergy distribution /(e) (given by the Fermi-Dirac distribu- 
tion in equilibrium) . Both are encoded in the interacting 
Green's functions through the relations 



Im G r (uj) = -it p(uj) , 

g K {u) = [2/H-i]i m g fl ( w ) 



(3) 



where G R = G++ - Q+~ and Q K = i[g++ + ]/2 are 
respectively the retarded and the Keldysh Green's func- 
tions. They obey the following Schwinger-Dyson equa- 
tions of motion 



g R (u) = [g*(u;)-i-z*(u,)]- 
g K (oj) = \g R (tj)\ 



(4) 



where T,^ K are the self-energy kernels corresponding to 
the Coulombic interaction. In order for the impurity 
problem to locally describe the lattice problem, one im- 
poses that it has the same (local) density of states and 
distribution function. This corresponds to the following 
self-consistent equations for the interacting Green's func- 
tions 

' g*(w) = G R (u;) , 
G K (u) = G K {u), 



(5) 
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FIG. 1: The impurity model consists in an impurity coupled to a 
set of leads with chemical potential Hi = I qEa, I 6 Z. In the frame- 
work of NESS-DMFT, the electric-field-driven Hubbard model in 
contact with a thermostat is solved through a self-consistent map- 
ping to the impurity model. 



where G R / K (u) = (a/2ir) 2 J dk G r / k (uj, k) are the local 
Green's functions and G R / K (uj, It) are the Green's func- 
tions of the lattice model. In our electric-field-driven 
setup, they correspond to the gauge-invariant Green's 
functions defined as G r / k (vo,k) in [8]. They are the 
solutions of the steady-state Schwinger-Dyson equations 
on the lattice, see Eqs. (2) and (3) in [1]. 

We solve for all of the lattice and impurity Green's 
functions with the NESS-DMFT algorithm described 
in Q. The self-energy kernels Yi^ K are now computed 
by solving the following impurity model with the steady- 
state impurity solver that we describe at the end of this 
Letter. The IPT solution is used as an initial guess. 

Impurity model. Let us re-interpret the Weiss fields 
as the result of integrating over the degrees of freedom 

of an effective environment described by the hybridiza- 

/? / a' 

tion kernels S C nv ■ The out-of-equilibrium nature of the 
environment is encoded in its energy distribution func- 
tion / cnv (e) defined by Ef nv = (2/ cnv - 1) Im E« nv . The 
non-interacting Green's functions of the impurity obey 
the following Schwinger-Dyson equations (assuming that 
a steady-state can be reached) 



g R {u) =[ w -s«»] 1 , 



(6) 



Similarly to the equilibrium case, there is no unique 
way of modeling the impurity environment since different 
environments can yield the same hybridization kernels. 
Nonetheless, a model motivated on solid physical ground 
will facilitate suitable solution schemes. Below, we first 
describe a generic impurity model, then we specialize it 
to the case of the electric-field-driven Hubbard model. 

Let us decompose the non-equilibrium environment 
into a collection of equilibrium non-interacting fcrmionic 
reservoirs, labeled by I, with different density of states 
pi (e) , temperatures 7] and chemical potentials pi . These 
reservoirs are linearly coupled to the impurity with the 
coupling constants i;. The corresponding hybridization 
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kernels read 

f Im S fnvM 



2T, 



(7) 



This decomposition is fairly generic and can be adapted 
to a wide class of steady-state impurity problems since 
it can accommodate any E^ lv (w) and T,^ nv (uj). We now 
specialize this approach to the case of the electric-field- 
driven Hubbard model. 

Let us consider for a while the potential gauge in which 
the electric field is rendered by a linear ramp of the scalar 
potential throughout the lattice: <fr(x) = — qE ■ x. If one 
singles out a site, its direct environment is composed of 
four neighbors: two are on the equi-potential (along the 
y direction) and the two others are shifted by ±qEa. All 
sites being equivalent, each of these neighbors feels the 
same kind of environment, and so on and so forth. 

Based on this discussion, we model the impurity en- 
vironment by a discrete "ladder" of equilibrium leads 
shifted in energy by Vi = I <p with tp = qEa and I G Z. 
A schematic representation of the model is depicted in 
Fig. [TJ We consider that the leads are at the same tem- 
perature Tl and have the same density of states pl (mod- 
ulo the energy shift) so that T\ = Tl, pi(e) = pl{c ~ Vi) 
and pi = Vi. In our symmetric case, we also have U = 
and we expect io/2 > t\ > t% > ... . This particular en- 
vironment yields a steady-state energy distribution func- 
tion 



/onv(e) 



(8) 



where /t l (e) = [1 + cxp(e/T£,)] _1 is the equilibrium 
Fermi-Dirac distribution at temperature Tl- Notice that 
although the distribution function / env is clearly not ther- 
mal for finite electric fields [see for instance Fig- EJb)] , Tl 
sets a unique thermal fluctuation scale in the problem. 
Below we validate the model by determining numerically 
the parameters of the model (the ti's, Tl and pl) that 
correctly reproduce the hybridization kernels (w). 

Noteworthy, the charge and energy local conservation 
laws read differently in the impurity problem than in the 
original lattice model. In the lattice, electrons flow along 
the x-axis. There is no net current of electrons between 
the lattice sites and the thermostat. However, there is 
finite heat current since the thermostat absorbs and dis- 
sipates the work that is performed by the electric field. 
On the impurity side, there are both an electronic cur- 
rent and a heat current between the impurity and each 
lead, and the dissipation takes place in every lead. 

This impurity model is strictly equivalent to an im- 
purity coupled to identical equilibrium reservoirs at the 
same chemical potential [i.e. pi = 0, pi(e) = Pl{^)] via 
some time-dependent couplings tie lVlt . This equivalence 
corresponds on the lattice side to the U(l) gauge invari- 
ance which states that the electric field can be rendered 
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FIG. 2: Fit of the environment: joint determination of PL(t), 
T L , t , and ti for U = 3 and E = 3 (T = 0.05, T = 0.3). (a) 
ImS^je) (plain line) is fitted with the dashed line which is the 
sum of the terms I = — 1 (red), 1 = (blue), and I = 1 (green) 
in Eq. Q. (b) / en v(e) (plain line) is fitted with the dashed line 
by the expression JS). High-energy features, |e| > \q\Ea, cannot 
be properly captured because of the truncation to the first three 
relevant leads. The other parameters of the fit are Tl ~ 0.05, 
t w 0.32 and t\ ks 0.18. 



by a linear ramp potential (potential gauge) or by some 
time-dependent hopping parameters (so-called temporal 
gauge). 

In practice, it is necessary to truncate the model to a 
finite number of leads: I = 0, ±1, ±2, .... ±l max . This ap- 
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proximation is valid as long as ^2=0* ^1 ^ Si^z 
Practically, we determine l m ax as the smallest integer 
such that /env(e) can be correctly fitted by Eq. ©. For 
simplicity, we choose to work with parameters E and U 
such that we can keep the leads I = —1,0,1 and dis- 
card the other leads (l m ax = 1)- This truncated im- 
purity model corresponds to an impurity coupled sym- 
metrically to two leads (/ = ±1) creating a voltage bias 
V1 — V-1 = 2qEa, and to one lead (/ = 0) at zero chemical 
potential po — Vq = 0. 

We determine the parameters of the model, to, t\, Tl, 
and yOi(e), by fitting the expressions of the hybridiza- 
tion kernels given in Eqs. (J7]) via a Newton's steepest de- 
scent method. To demonstrate the validity of the model, 
we give in Fig. [2] the result of the fitting procedure for 
U = 3 and E = 3 after the full DMFT solution has 
converged (see below). In Fig. [5Ja), we plot the decom- 
position of Im S^ w (e) as the sum of t 2 pL(e — tp), t^pL^) 
and t\pL{t + p)- The local maxima of Im T, env (e) at 
e = ±qEa and e = naturally come from the replication 
of unique maximum of pl(^) at e = 0. In Fig. EJb), we 
plot the corresponding energy distribution function. The 
jump around e = is accounted by the / = lead while 
the two symmetric jumps around e = ±qEa are described 
by the leads I = ±1. As a result of the truncation, the 
high-energy features above |e| > \q\Ea are not properly 
captured and t\ is slightly overshot. 

Impurity solver. To demonstrate the practical rele- 
vance of the impurity model we introduced, we solve it 
and compute non-perturbatively the interaction contri- 
bution to the retarded self-energies by generalizing 
to multiple leads the steady-state impurity solver that 
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self-energy is obtained by 



FIG. 3: (a) Retarded self-energy Im T,§(lo) obtained by the NESS- 
HF solver, (b) Local density of states p(e). The dashed lines cor- 
respond to the IPT solution. (U = 3, E = 3, T = 0.05, T = 0.3). 



was recently developed by Han and Heary in the context 
of a two- lead environment pLUlj. Unlike the IPT solver, 
this solver provides ^-derivable solutions and is therefore 
suited to cases away from half-filling. 

Building on Hershfield's expression for the steady-state 
density matrix [12| , Han and Heary gave an effective Mat- 
subara description of the steady state. This allows the 
use of standard many-body equilibrium tools to tackle 
the strong interaction at the cost of introducing some 
imaginary chemical potentials: itp m = i2m7rTi, m G 7L. 
After the impurity problem is solved in imaginary time 
for each m, the solutions are analytically continued to 
the real-time real- voltage domain. 

We refer the reader to Q for the technical details as 
we follow closely the steps explained there. The non- 
interacting impurity Green's function in Matsubara fre- 
quency and imaginary-voltage reads: 



z=o,±i 



Ep^Pi( £ ~ V p) Z nm.l ~ 



(9) 



with z nm ,i = \uj n - l(iip m - <p), iuj n = i(2n + 1)ttT l , and 
p (e) = -Im G^(e)/w. The differences with Eq. (32) 
in [5[ lie in the presence of the I = term and the 
frequency-dependent environment which is determined 
self-consistently. 

The retarded self-energies S^(io-' 7] ,. U£ m ) are obtained 
by means of a Hirsch-Fye algorithm [13j . More sophisti- 
cated algorithms have already been used in the context 
of a two- lead static environment [13]. S^(w) is then ob- 
tained after a double analytical continuation: iw„ — > u) 
and iip m — > (p. This is performed numerically by fitting 
all the Tj^j (iw„, iip m ) as 



£y(iu; n ,i^ m ) 



E 

7£Z 



dc 



cr 7 (e)Q 7 (e,i^ m - tp) 



(10) 



where we truncate the sum at the order 7 max = 5 and 
the function Q 7 (e, z) is fitted by the simple Pade ap- 
proximant ^^ T |^ . The parameters of the fit are the 
functions er 7 , C 7 and D 1 . 



Im Ey(w) = — 7T G-fi^) ■ 



(11) 



We used a simple minimization procedure based on New- 
ton's steepest descent method. Although the success of 
analytical continuations is never exempt from some nu- 
merical 'cooking' techniques, this double analytical con- 
tinuation was unexpectedly rather easy to perform, per- 
haps thanks to the presence of some finite dissipation 
in the bulk that moves dangerous poles away from the 
real axis. In Fig. [3Ja) we plot Im S^(cj) for U = 3 and 
E = 3 and compare with the result obtained with the 
IPT solver. The sharper edges around u> = are respon- 
sible for a better defined quasi-particle peak in the local 
density of states [see Fig. EJb)] but altogether, the re- 
sults are quite close to the ones obtained with IPT. This 
agreement for relatively small U validates our impurity 
model and impurity solver. 

We approximate the steady-state distribution function 
of the interacting impurity by the one of its environment: 
/(e) = /cnv(e)- It is a reasonable approximation for the 
following reasons: (i) it is exact in equilibrium, (ii) it 
is exact in the non-interacting case, (iii) it is consistent 
with charge and energy conservations, In turn, this gives 
us a simple way to estimate the Keldysh component of 
the self-energy: 



E£(w) = [2/(w)-l]Im E£(w) 



(12) 



Once the impurity is solved, the self-energy kernels Y>^ K 
are used to recompute the lattice Green's functions, a 
new impurity problem is defined and solved again, so on 
and so forth until convergence is reached. 

Conclusion. Although the work presented here is cen- 
tered around the electric-ficld-drivcn Hubbard model, the 
impurity model, the dictionary between the lattice side 
and the impurity side, and the steady-state impurity 
solver allow to consider further questions, such as the 
effect of chemical substitution (doping) and pressure ef- 
fects under an electric field, not achievable with simpler 
perturbation methods. Also, not only the steady-state 
impurity solver enables to study the transport properties 
of quantum dots driven by leads with a realistic density 
of states, it can also address fundamental questions such 
as the effect of a quantum critical point (gapped leads) 
on the out-of-equilibrium Kondo physics |15| . 
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The real-frequency retarded 
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